Departamento de Ecologia IB-USP

Um problema

Qual é a alturas média das årvores em uma floresta?

InferĂȘncia estatĂ­stica

  • Fazer afirmaçÔes sobre as caracterĂ­sticas de uma população, com base em informaçÔes dadas por amostras.

  • População Ă© qualquer conjunto de elementos que vocĂȘ estĂĄ observando

  • Amostra Ă© qualquer subconjunto desta população.

Ajustando um modelo

  • Estimar parĂąmetros da população baseados na amostra, modelo teĂłrico de distribuição de probabilidades

Modelo de distribuição para altura das årvores

  • VariĂĄvel AleatĂłria de distribuição de probabilidades Normal com mĂ©dia (ÎŒ) e desvio padrĂŁo (σ).

## [1] 22.72789
## [1] 32.3065
## [1] 5.683881

Definindo variĂĄveis aleatĂłrias

A variĂĄvel aleatĂłria X Ă© uma variĂĄvel que tem um valor Ășnico (determinado aleatoriamente) para cada resultado de um experimento (lembre-se esse termo veio da teoria de probabilidade).

Exemplos:

  1. Altura das ĂĄrvores de uma floresta (!!)

  2. NĂșmero de presas capturadas por um predador em um determinado dia;

  3. Comprimento de um peixe adulto selecionado aleatoriamente.

As variĂĄveis aleatĂłrias podem ser discretas ou contĂ­nuas.

VariĂĄvel aleatĂłria discreta

NĂșmero ou a quantidade observada na unidade experimental ou tentativa.

  • Representada por nĂșmeros inteiros (0, 1, 2, 3, 4…);

  • NĂŁo pode conter nĂșmeros negativos;

  • NĂșmero finito de possibilidades;

  • Podemos achar a probabilidade de cada evento.

VariĂĄvel aleatĂłria contĂ­nua

Ex: peso, altura, distĂąncia, pH, biomassa, etc.

  • Representada por nĂșmeros nĂŁo inteiros (1.3, - 1.54, - 1.7);

  • Pode conter nĂșmeros negativos;

  • NĂșmero infinito de possibilidades;

  • Probabilidade de cada evento Ă© zero.

Função de probabilidade

Associa cada possĂ­vel valor da variĂĄvel (X) Ă  sua probabilidade de ocorrĂȘncia P(X).

Quando conhecemos todos os valores de uma variåvel aleatória, juntamente com suas respectivas probabilidades, temos uma distribuição de probabilidades.

Algumas DistribuiçÔes de Probabilidade

Bernouli: X ~ Bernouli(p)

  • Jogar moedas: cara 0, coroa 1

  • Presença/ ausĂȘncia de uma espĂ©cie em locais amostrados: 0 - 1

  • ParĂąmetro: p - probabilidade de sucesso

E[X] = p

var[X] = p(1-p)

Binomial: X ~ Bin(n, p)

  • Distrib. Discreta do nĂșmero de sucessos em uma sequĂȘncia de n tentativas

E[X] = np

e a variĂąncia Ă©

var[X] = np(1-p)

  • NĂșmero sementes predadas em cada experimento contendo 20 semente.

  • NĂșmero de animais infectados em relação ao nĂșmero de capturados em cada local

Exemplo Binomial

  • Probabilidade de predação de girinos por odonatas. P = 0.30

  • Temos 6 girinos num lago.

  • Qual robabilidades de que 0, 1, 2, 3, 5 ou todos sejam predados?

Poisson: X ~ Pois(λ)

  • Distribuição de probabilidade discreta.

  • Probabilidade de uma sĂ©rie de eventos ocorrem em um perĂ­odo fixo de tempo, ĂĄrea, volume, quadrante, etc.

E[X] = (λ)

var[X] = (λ)

  • AbundĂąncia de espĂ©cies em um local

  • NĂșmero de indivĂ­duos capturados por tempo

Exemplo Poisson

  • NĂșmero de visitas de borboletas em uma planta em 15 min: 10
  • Probabi A probabilidade de uma borboleta visitar Ă© a mesma para quaisquer dois perĂ­odos de tempo de igual comprimento. Abaixo o histograma dessa distribuição de probabilidade.

Normal: X ~ N(ÎŒ, σ)

  • SimĂ©trica, contĂ­nua

– Metade (50%) estĂĄ acima (e abaixo) da mĂ©dia

– Aproximadamente 68% estĂĄ dentro de 1 desvio padrĂŁo da mĂ©dia

– Aproximadamente 95% estĂĄ dentro de 2 desvios padrĂ”es da mĂ©dia

– Virtualmente todos os valores estĂŁo dentro de 3 desvios padrĂ”es da mĂ©dia

Teorema Central do Limite

"toda soma de variĂĄveis aleatĂłrias independentes de mĂ©dia finita e variĂąncia limitada Ă© aproximadamente Normal, desde que o nĂșmero de termos da soma seja suficientemente grande".

Independentemente do tipo de distribuição da população, na medida em que o tamanho da amostra aumenta, a distribuição das médias amostrais tende a uma distribuição Normal.

Exemplo Normal

  • Qual Ă© a probabilidade de que um peixe capturado aleatoriamente tenha 20,15 cm ou mais?

  • MĂ©dia da população Ă© 17,1 cm e o desvio padrĂŁo Ă© de 1,21 cm?

prob <- pnorm(q = 20.15, mean = 17.1, sd = 1.21, lower.tail = F)
prob
## [1] 0.005856729

Voltando ao nosso problema

  • Ddados de duas espĂ©cies
levels(arv$sp)
## [1] "sp1" "sp2"

Serå que a altura média de uma espécie nessa floresta é diferente da altura média da outra espécie?

Ou seja, serå que as espécies pertencem a uma mesma distribuição e as diferenças encontradas são devido ao acaso ou serå que cada espécie é uma variåvel aleatória com médias diferentes?

DistribuiçÔes de frequĂȘncias das alturas para cada espĂ©cie

library(lattice)
histogram(~alt|sp, 
          data = arv,
          layout = c(1,2),
          strip=F,
          strip.left=T)

Diferença entre as médias das alturas

mean(arv$alt[arv$sp=="sp2"]) - mean(arv$alt[arv$sp=="sp1"])
## [1] 5.109466

Posso afirmar que as médias são diferentes?

HipĂłteses sobre as alturas

H0 - não hå diferença entre as alturas das duas espécies. A média da altura da espécie 1 é igual à média da altura da espécie 2 e a diferença observada se deve ao acaso.

H1 - hå diferença entre as alturas das duas espécies. A média da altura da espécie 1 é diferente da média da altura da espécie 2.

Erros associados a escolher cada hipĂłtese como verdadeira

AFIRMAÇÃO: nĂŁo hĂĄ diferença entre as alturas das espĂ©cies,
VERDADE: hå diferenças.

AFIRMAÇÃO: hĂĄ diferença entre as alturas das espĂ©cies,
VERDADE: não hå diferenças.

Definição formal dos erros:

ERRO DO TIPO I (α): rejeitar a hipótese nula dada que ela é verdadeira.

ERRO DO TIPO II: aceitar a hipĂłtese nula dado que ela Ă© falsa.

Testando as hipĂłteses

Contrastar a probabilidade de que a diferença que vocĂȘ encontrou nas mĂ©dias das alturas das duas espĂ©cies pode se considerada devido ao acaso ou nĂŁo.

Qual é a probabilidade de que as duas amostras pertencem à mesma população de medidas?

Podemos ou não rejeitar a hipótese nula (de que a diferença é ao acaso), baseado em alguma probabilidade de estarmos cometendo algum erro (α).

Graus de liberdade

"O nĂșmero de observaçÔes independentes menos o nĂșmero de parĂąmetros estimados."

Quantas observaçÔes independentes foram usadas para calcular a média e a variùncia dos dados sob a hipótese nula?

100 ĂĄrvores
2 parĂąmetros
98 graus de liberdade

Teste t

Não entraremos em detalhes de como se calcula a estatística t e como é a sua distribuição de probabilidades (recomendo ver na wikipedia). O importante a saber aqui é que tendo as duas amostras modeladas como uma distribuição normal, calculamos a estatística t e com um certo valor de graus de liberdade, nós podemos recorrer às tabelas da distribuição t para ver qual a probabilidade de se obter aquele valor que obtivemos dada que a hipótese nula é verdadeira. Se esta probabilidade for muito pequena, a chance de estarmos caindo no erro do tipo 1 é pequena.

No R, a função usada para fazer um teste t é t.test.

t.test(arv$alt~arv$sp, var.equal=T)
## 
##  Two Sample t-test
## 
## data:  arv$alt by arv$sp
## t = -5.0125, df = 98, p-value = 2.387e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.132311 -3.086621
## sample estimates:
## mean in group sp1 mean in group sp2 
##          20.17316          25.28262
#o argumento var.equa=T diz que as variùncias das duas espécies são as mesmas
# podemos fazer o teste com a premissa de que as variĂąncias nĂŁo sejam iguais,
# usando var.equal=F, isso vai "consumir" graus de liberdade.

Vamos interpretar este resultado: temos o valor de t calculado e os graus de liberade (jå calculado pela gente antes). Com estes 2 valores, o que o R faz é calcular a probabilidade de que a hipótese nula seja verdadeira. Esse é o valor-p, que neste exemplo deu um valor beem pequeno. Vemos então que podemos rejeitar H0, com uma probabilidade muito baixa de estarmos incorrendo no erro de tipo 1, e desta forma acreditamos que as espécies de årvores tenham distribuiçÔes de alturas diferentes (médias diferentes).

Modelos lineares

E se tivéssemos na verdade 3 espécies de årvores para compararmos as alturas? Acabamos de ver que o teste t serve apenas para comparar duas amostras. O que chamamos de modelos lineares é uma grande família de modelos estatísticos que modelam uma relação linear entre a variåvel dependente (Y) e a(s) variåvel(is) independente(s). Nesta família encontra-se o teste t que acabamos de ver, a Anålise de Variùnci (ANOVA), a regressão, etc.

Chamamos de Modelos Lineare Gerais (LMs), aqueles que possuem como premissa a variåvel dependente Y como de distribuição Normal, e de Modelos Lineares Genearlizados (GLMs) os que assumem outros tipos de distribuição para a variåvel dependente (como a Binomial e Poisson). No R, a função que usamos para estes modelos são lm e glm (hå alguns mais específicos como o prórpio t.test e o aov para ANOVA). O R lida com todos os modelos lineares de maneira similar, e um detalhe importante é que todas as anålises desta família devem ser salvas em um objeto, para seu resultado ser apresentado após o comando summary.

É importante lembrar que neste roteiro estamos lidando com modelos paramĂ©tricos, ou seja, que modela a variĂĄvel de interesse como uma variĂĄvel aleatĂłria pertencente a uma certa distribuição de probabilidades. Existem outras abordagens estatĂ­sticas nĂŁo-paramĂ©tricas que fazem testes que nĂŁo tem como premissa a distribuição de probabilidades. Veja um exemplo de teste t nĂŁo paramĂ©trico neste roteiro.

Mesmo fazendo parte de uma mesma família de modelos, nas sessÔes seguintes vamos separar os modelos lineares pelos nomes tradicionais das anålises: anålise de variùncia, regressão linear simples e anålise de covariùncia. Depois, veremos alguns exemplos de GLMs aplicåveis a dados ecológicos.

ANOVA

Vamos entĂŁo incluir no nosso exemplo inicial uma terceira espĂ©cie amostrada (e vamos diminuir o nĂșmero de amostras por espĂ©cie para ficar mais fĂĄcil de fazer os cĂĄlculos passo a passo):

alt.sp3 <- c(36.90076761,33.23131727,32.82767311,24.93410577,31.87626329,
         28.76824248,26.6436144,31.41238398,27.65383929,31.70425719,
         33.60752445,28.30263354,29.32210674,23.4790054,22.20793046,
         18.28797293,23.39044736,28.75820917,29.94344703,34.47430235,
         32.09820862,26.99845592,27.64858951,31.90557306,33.63236368,
         25.35966179,35.00885172,32.89125598,18.98802229,30.8516906,
         22.04191398,35.37232889,35.67150408,34.30257037,24.43038475,
         20.6274663,31.15717916,23.27023231,36.20568545,34.33068448,
         28.40769054,31.00578874,29.50989158,31.69915991,37.25975797,
         30.1279997,38.623371,35.36472829,29.88304259,31.00145491)

sp3 <- data.frame(alt=alt.sp3, sp=rep("sp3",50))
arv2 <- rbind(arv,sp3)

amostras <- seq(1,150,by=5)

arv3 <- arv2[ amostras, ]

Neste exemplo lidaremos com um modelo simples de anova chamado one way ANOVA, porque lida apenas com uma variĂĄvel independente categĂłrica (fator).

A ANOVA testarå a hipótese nula de que as médias das alturas das 3 espécies não diferem.

Faremos agora a ANOVA "na unha", para entendermos melhor cada passo da anålise e a interpretação de seus resultados para depois usarmos a função do R que faz esta anålise diretamente.

ANOVA na unha

Vamos calcular os valores da tabela de ANOVA. Começando com os desvios quadråticos, ou seja, quanto os dados desviam da média (idéia parecida com a variùncia). O ponto importante é que essa variação é aditiva e portanto, pode se decomposta. A variação total é decomposta em variação relacionada ao tratamento (entre grupos), no nosso caso às espécies, e uma variação interna dos grupos (chamada de erro). A estatística F é a razão entre essas variaÔes após dividir cada uma delas pelos seus respectivos graus de liberdade. Complicou? Vamos fazer os cålculos e ver os gråficos para ver se entendemos melhor:

A tabela de ANOVA que vamos construir Ă© essa:

Primeiro vamos mudar a forma do nosso data frame para facilitar os comandos

arv4 <- data.frame(sp1=arv3$alt[arv3$sp=="sp1"],
                   sp2=arv3$alt[arv3$sp=="sp2"],
                   sp3=arv3$alt[arv3$sp=="sp3"])
arv4
##          sp1      sp2      sp3
## 1  21.282000 19.14508 36.90077
## 2  19.875790 18.66252 28.76824
## 3  23.562570 25.10581 33.60752
## 4  19.313744 18.10266 18.28797
## 5  26.308616 21.48084 32.09821
## 6   8.301795 28.08559 25.35966
## 7  27.192502 20.47805 22.04191
## 8  15.941866 29.22927 20.62747
## 9  20.596624 19.32911 28.40769
## 10 22.048945 18.21446 30.12800
boxplot(arv4)

Vamos calcular a média geral, as médias para cada espécie, e as diferenças entre a média geral e a média para cada espécie, para chegar ao valor da soma dos quadrados total:

media.geral <- mean(arv3$alt)
media.geral
## [1] 23.28284
medias.sps <- apply(arv4, 2, mean)
medias.sps
##      sp1      sp2      sp3 
## 20.44245 21.78334 27.62274
dif.geral <- arv4 - media.geral
ss.especies <- dif.geral^2
ss.especies
##                sp1       sp2        sp3
##  [1,]   4.00337202 17.121094 185.447876
##  [2,]  11.60801175 21.347420  30.089610
##  [3,]   0.07824755  3.323227 106.599051
##  [4,]  15.75374565 26.834280  24.948725
##  [5,]   9.15530094  3.247214  77.710675
##  [6,] 224.43179074 23.066347   4.313177
##  [7,]  15.28543337  7.866884   1.539904
##  [8,]  53.88994194 35.359990   7.051024
##  [9,]   7.21577099 15.631997  26.264064
## [10,]   1.52250306 25.688498  46.856173
ss.total <- sum(ss.especies)
ss.total
## [1] 1033.251

Para calcular os desvios quadråticos totais, nós subtraímos cada altura das årvores pela média geral e elevamos ao quadrado. Veja o gråfico:

vetor.cor <- rep(1:3, each=10) #vetor de cores

plot(x = c(1:30), y = arv3$alt, ylim=c(10,40),
     pch=(rep(c(15,16,17), each=10)),
     col=vetor.cor,
     ylab="Variåvel Resposta", xlab="ObservaçÔes",
     main="Variação total")
    for(i in 1:30)
    {
    lines(c(i,i),c(arv3$alt[i],mean(arv3$alt)),col=vetor.cor[i])
    }
    abline(h=media.geral)

Agora vamos fazer a somatĂłria dos desvios quadrĂĄticos dentro de cada grupo (ss.intra).

O grĂĄfico para entender esse cĂĄlculo:

vetor.medias<-rep(medias.sps, each=10)

plot(c(1:30), arv3$alt, ylim=c(10,40),
     pch=(rep(c(15,16,17),each=10)),
     col=vetor.cor,
     main="Variação Intra Grupos",
     ylab="Variåvel Resposta", xlab="ObservaçÔes")
    for(i in 1:30)
    {
    lines(c(i,i),c(vetor.medias[i],arv3$alt[i]),col=vetor.cor[i])
    }
    lines(c(1,10),c(medias.sps[1],medias.sps[1]),col=1)
    lines(c(11,20),c(medias.sps[2],medias.sps[2]),col=2)
    lines(c(21,30),c(medias.sps[3],medias.sps[3]),col=3)

CĂĄlculo dos valores:

medias.sps
##      sp1      sp2      sp3 
## 20.44245 21.78334 27.62274
ss.sp1=sum((arv4$sp1-medias.sps["sp1"])^2)
ss.sp1
## [1] 262.2655
ss.sp2=sum((arv4$sp2-medias.sps["sp2"])^2)
ss.sp2
## [1] 157.0018
ss.sp3=sum((arv4$sp3-medias.sps["sp3"])^2)
ss.sp3
## [1] 322.4728
ss.intra=ss.sp1+ss.sp2+ss.sp3
ss.intra
## [1] 741.7401

A soma dos quadrados entre grupos:

plot(c(1:30), vetor.medias, ylim=c(10,40), 
  pch=(rep(c(15,16,17),each=10)),
  col=vetor.cor, 
  main="Variação Entre Grupos", 
  ylab="Variåvel Resposta", xlab="ObservaçÔes")
 for(i in 1:30)
    {
    lines(c(i,i),c(vetor.medias[i],mean(vetor.medias)),col=vetor.cor[i])
    }
    abline(h=media.geral)
    points(c(1:30),arv3$alt, ylim=c(10,50), 
           pch=(rep(c(0,1,2),each=10)), col=vetor.cor, cex=0.5)

medias.sps
##      sp1      sp2      sp3 
## 20.44245 21.78334 27.62274
media.geral
## [1] 23.28284
ss.entre=10*sum((medias.sps-media.geral)^2)
ss.entre
## [1] 291.5112
#conferindo os cĂĄlculos
ss.intra+ss.entre
## [1] 1033.251
ss.total
## [1] 1033.251

CĂĄlculo do F

# Desvios médios
ms.entre=ss.entre/2
ms.intra=ss.intra/27
ms.entre
## [1] 145.7556
ms.intra
## [1] 27.47186
#F - razĂŁo das variĂąncias
F.sps=ms.entre/ms.intra
F.sps
## [1] 5.305634

Distribuição F

Vamos ver na distribuição F qual a probabilidade de termos econtrado o valor F.sps ao acaso:

curve(expr=df(x, 2,27),
      main="Distribuição F de Fisher (df=2,27)", 
      xlab="Valor F",
      ylab="Densidade ProbabilĂ­stica (df)",
      xlim=c(2,12))
abline(v=F.sps, col="red")
abline(h=0, lty=2)

xf=seq(F.sps, 12, 0.01)
ydf=df(xf, 2, 27)
polygon(c(F.sps,xf),c(0,ydf),col="red")

text(x= 7,y=0.08,paste("pf(x) =",
     round(pf(F.sps,2,27,lower.tail=F),4)), 
     cex=0.8, col="red")

CĂĄlculo do P

p.sps=pf(F.sps, 2, 27, lower.tail=FALSE)
p.sps
## [1] 0.01139248

Para mais informaçÔes sobre o F e as comparaçÔes com o teste t, veja esse roteiro.

A tabela final

Colocando os dados calculados em nossa tabela de anova

ANOVA no R

A função que constrói esta tabela de ANOVA e faz o teste estatístico é a aov. Vamos comparar:

anova.sps <- aov(alt~sp, data=arv3)
summary(anova.sps)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## sp           2  291.5  145.76   5.306 0.0114 *
## Residuals   27  741.7   27.47                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Agora nĂłs sabemos de ondem vem os nĂșmeros nesta tabela e como foi calculado o P, certo?

Testes Ă  posteriori

Onde estão as diferenças?

  • veja os grĂĄficos

  • faça teste Ă  posteriori - Tukey HDS

TukeyHSD(anova.sps)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = alt ~ sp, data = arv3)
## 
## $sp
##             diff        lwr       upr     p adj
## sp2-sp1 1.340893 -4.4708805  7.152667 0.8360276
## sp3-sp1 7.180300  1.3685259 12.992073 0.0132022
## sp3-sp2 5.839406  0.0276327 11.651180 0.0487468
#revendo grĂĄfico dos dados
stripchart( arv3$alt~arv3$sp,
           vertical = TRUE, pch = 16, 
           col = c("black", "red", "green"))

#boxplot
boxplot(arv4)

Outra maneira de fazer a ANOVA no R:

lm.anova.sp <- lm(alt~sp, data=arv3)
summary.aov(lm.anova.sp) #mesmo que o summary do aov
##             Df Sum Sq Mean Sq F value Pr(>F)  
## sp           2  291.5  145.76   5.306 0.0114 *
## Residuals   27  741.7   27.47                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mostra o resultado de forma um pouco diferente, consegue entender?
summary(lm.anova.sp) 
## 
## Call:
## lm(formula = alt ~ sp, data = arv3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.1407  -3.0002  -0.0742   3.2719   9.2780 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   20.442      1.657  12.334 1.32e-12 ***
## spsp2          1.341      2.344   0.572  0.57202    
## spsp3          7.180      2.344   3.063  0.00492 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.241 on 27 degrees of freedom
## Multiple R-squared:  0.2821, Adjusted R-squared:  0.229 
## F-statistic: 5.306 on 2 and 27 DF,  p-value: 0.01139

RegressĂŁo linear simples

  • Y contĂ­nuo
  • X continuo
plot(arv3$alt ~ nutri)

Modelo da RegressĂŁo

alt.nutr <- lm(alt ~ nutri, data = arv3)
summary(alt.nutr)
## 
## Call:
## lm(formula = alt ~ nutri, data = arv3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7176 -1.7171  0.0088  1.4124  9.2072 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 14.38711    1.05832  13.594 7.44e-14 ***
## nutri        0.58006    0.05971   9.714 1.82e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.906 on 28 degrees of freedom
## Multiple R-squared:  0.7712, Adjusted R-squared:  0.763 
## F-statistic: 94.37 on 1 and 28 DF,  p-value: 1.818e-10

Plot da RegressĂŁo

plot(arv3$alt ~ nutri)
abline(alt.nutr, col = "red") # reta da regressĂŁo

#reta da hipĂłtese nula de ausĂȘncia de efeito dos nutrientes
abline(h = mean(arv3$alt), col = "blue", lty = 2)  

AnĂĄlise de CovariĂąncia

  • Y continuo

  • X1 contĂ­nuo

  • X2 categĂłrico

plot(alt~nutri, data=arv3, col=arv3$sp)

Modelo da ANCOVA: aditivo

alt.sp.nutr <- lm(alt ~ nutri + sp, data = arv3)
summary(alt.sp.nutr)
## 
## Call:
## lm(formula = alt ~ nutri + sp, data = arv3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6682 -0.9554 -0.1585  1.3581  7.6134 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  13.3971     1.1034  12.142 3.23e-12 ***
## nutri         0.5256     0.0561   9.369 8.09e-10 ***
## spsp2         1.6421     1.1423   1.437  0.16251    
## spsp3         3.8326     1.1965   3.203  0.00357 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.553 on 26 degrees of freedom
## Multiple R-squared:  0.836,  Adjusted R-squared:  0.817 
## F-statistic: 44.16 on 3 and 26 DF,  p-value: 2.401e-10

GrĂĄfico ANCOVA aditivo

Modelo ANCOVA com interação

alt.sp.nutr2 <- lm(alt ~ nutri + sp + nutri:sp, data=arv3)

alt.sp.nutr3 <- lm(alt ~ nutri*sp, data=arv3)
summary(alt.sp.nutr3)
## 
## Call:
## lm(formula = alt ~ nutri * sp, data = arv3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6724 -0.9444 -0.3313  1.0416  7.0702 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13.40167    1.42727   9.390 1.66e-09 ***
## nutri        0.52527    0.08906   5.898 4.38e-06 ***
## spsp2        2.93717    1.97080   1.490    0.149    
## spsp3        0.43636    2.75666   0.158    0.876    
## nutri:spsp2 -0.10095    0.12423  -0.813    0.424    
## nutri:spsp3  0.17187    0.14350   1.198    0.243    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.474 on 24 degrees of freedom
## Multiple R-squared:  0.8578, Adjusted R-squared:  0.8282 
## F-statistic: 28.96 on 5 and 24 DF,  p-value: 2.003e-09

Gråfico ANCOVA com interação

Tabela de ANOVA do modelo

summary.aov(alt.sp.nutr3)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## nutri        1  796.8   796.8 130.179 3.53e-11 ***
## sp           2   66.9    33.5   5.467   0.0111 *  
## nutri:sp     2   22.6    11.3   1.846   0.1796    
## Residuals   24  146.9     6.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelos lineares generalizados

  • Outras distribuiçÔes de probabilidades alĂ©m da Normal

  • VariĂąncias nĂŁo homogĂȘneas

Mais comuns:

  • Poisson – Ășteis para dados de contagem

  • Binomial – Ășteis para dados com proporçÔes, ou presença/ausĂȘncia

Etapas do GLM

  1. Uma hipótese sobre a distribuição da variåvel resposta Yi.

  2. Especificação do modelo (variåveis X).

3.Função de ligação entre a média Yi e as variåveis X. Função que lineariza a regressão

FunçÔes de ligação

Distribuição link function
normal identity
poisson log
binomial logit
Gamma reciprocal

GLM Poisson

  • Y - NĂșmero de atropelamentos de anuros em uma estrada

  • X - distĂąncia a um Parque Nacional.

plot( atrop$TOT.N ~ atrop$D.PARK,
      xlab = " DistĂąncia ao parque",
      ylab = "NĂșmero de atropelamentos")

Modelo GLM Poisson

mod <- glm(TOT.N ~ D.PARK, data = atrop, 
           family = poisson)

summary(mod)
## 
## Call:
## glm(formula = TOT.N ~ D.PARK, family = poisson, data = atrop)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -8.1100  -1.6950  -0.4708   1.4206   7.3337  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.316e+00  4.322e-02   99.87   <2e-16 ***
## D.PARK      -1.059e-04  4.387e-06  -24.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1071.4  on 51  degrees of freedom
## Residual deviance:  390.9  on 50  degrees of freedom
## AIC: 634.29
## 
## Number of Fisher Scoring iterations: 4

Testando contra hipĂłtese nula

anova(mod, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: TOT.N
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                      51     1071.4              
## D.PARK  1   680.55        50      390.9 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

GrĂĄfico Modelo GLM Poisson

##   (Intercept)        D.PARK 
##  4.3164849800 -0.0001058506

GLM Binomial

  • Y - Quantidade de veados infectados dentre os capturados

  • X - quantidade de vegetação aberta

plot(prop.veados ~ veados$openland,
     ylim=c(0,1),
     xlab = "Porcentagem de vegetação aberta",
     ylab="Proporção de veados infectados")

Modelo GLM binomial

nao.infectado <- veados$amostrado- veados$infectado
y <- cbind(veados$infectado, nao.infectado)
mod.veado <- glm(y ~  openland, data=veados, family="binomial")
summary(mod.veado)
## 
## Call:
## glm(formula = y ~ openland, family = "binomial", data = veados)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.8052  -1.1463   0.6707   1.3416   3.5483  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.4820     0.1736   14.29   <2e-16 ***
## openland     -6.0336     0.4432  -13.61   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 376.18  on 31  degrees of freedom
## Residual deviance: 121.30  on 30  degrees of freedom
## AIC: 204.65
## 
## Number of Fisher Scoring iterations: 5

Teste GLM binomial

anova(mod.veado, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: y
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                        31     376.18              
## openland  1   254.88        30     121.30 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

GrĂĄfico GLM binomial

## (Intercept)    openland 
##    2.482035   -6.033583

Checando as premissas dos modelos

Teste de normalidade

Teste de Shapiro-Wilks

  • H0: dados sĂŁo normais
shapiro.test(arv3$alt)
## 
##  Shapiro-Wilk normality test
## 
## data:  arv3$alt
## W = 0.96262, p-value = 0.3609
shapiro.test(arv$alt[arv$sp == "sp1"])
## 
##  Shapiro-Wilk normality test
## 
## data:  arv$alt[arv$sp == "sp1"]
## W = 0.9906, p-value = 0.9591
shapiro.test(arv$alt[arv$sp == "sp2"])
## 
##  Shapiro-Wilk normality test
## 
## data:  arv$alt[arv$sp == "sp2"]
## W = 0.97453, p-value = 0.3502

Teste de homogeneidade de variĂąncias

Teste de Bartlett:

bartlett.test(arv3$alt~arv3$sp)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  arv3$alt by arv3$sp
## Bartlett's K-squared = 1.1109, df = 2, p-value = 0.5738

Teste de Levene:

leveneTest(arv3$alt~arv3$sp)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2   0.512  0.605
##       27

Checagem do ajuste do modelo - ResĂ­duos

RESÍDUO = Y(observado) - Y(dados ajustados)

A premissa de normalidade na verdade recai sobre os resĂ­duos do modelo.

resid(anova.sps)
##           1           6          11          16          21          26 
##   0.8395548  -0.5666556   3.1201253  -1.1287012   5.8661704 -12.1406501 
##          31          36          41          46          51          56 
##   6.7500566  -4.5005793   0.1541789   1.6065001  -2.6382600  -3.1208224 
##          61          66          71          76          81          86 
##   3.3224765  -3.6806771  -0.3024985   6.3022481  -1.3052922   7.4459310 
##          91          96         101         106         111         116 
##  -2.4542276  -3.5688778   9.2780228   1.1454976   5.9847796  -9.3347719 
##         121         126         131         136         141         146 
##   4.4754638  -2.2630830  -5.5808309  -6.9952785   0.7849457   2.5052549
#residuals(anova.sps) # o mesmo que a função de cima

plot do modelo

ResĂ­duos regressĂŁo

# regressĂŁo: alt ~ nutri
plot(alt.nutr)

ResĂ­duos ANCOVA aditivo

# ANCOVA modelo aditivo: alt ~ nutri + sps
plot(alt.sp.nutr)

ResĂ­duos ANCOVA interativo

# ANCOVA modelo interativo: alt ~ nutri * sps
plot(alt.sp.nutr2)

ResĂ­duos GLM poisson

# GLM poisson: atrop ~ dist.park
plot(mod)

ResĂ­duos GLM binomial

# GLM binoimial: prop.infectado ~ openland
plot(mod.veado)